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ABSTRACT 


The  spherical  earth  diffraction  subroutine  SPH35  in  the 
radar  propagation  code  SEKE  has  been  known  to  cause  errors  in 
propagation  loss  computations  for  a  range  of  combinations  of 
antenna  and  target  heights.  In  this  report  an  efficient  method 
to  evaluate  the  Airy  function  in  the  complex  plane  is 
presented.  This  method  uses  the  power  series  expansion  near 
the  origin  and  an  integral  representation  elsewhere.  It  is 
more  accurate  and  as  fast  as  the  method  employed  in  the  spheri¬ 
cal  earth  diffraction  subroutine  SPH35  that  evaluates  every 
Airy  function  of  Fock' s  series  by  a  fourth-order  polynomial  fit 
to  its  logarithm.  The  algorithm  presented  was  incorporated  in 
a  new  spherical  earth  diffraction  subroutine  (SPH35N) .  It  was 
found  that,  if  SEKE  uses  this  subroutine,  no  problems  arise  for 
normalized  heights  of  up  to  5000  (i.e.  about  350  km  at  VHF  or 
17  km  at  Ku  band) . 

The  subroutine  SPH35N,  described  in  this  report,  has  been 
used  in  the  versions  of  SEKE  running  at  Lincoln  Laboratory,  and 
is  in  the  version  of  SEKE  currently  being  supplied  to  other 


users . 
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1.  INTRODUCTION 

SEKE  is  a  computer  model  that  predicts  the  one-way 
propagation  factor  over  irregular  terrain  by  selecting,  based 
on  terrain  geometry,  one  algorithm  or  a  weighted  combination  of 
algorithms  designed  to  compute  specular  reflection,  spherical 
earth  diffraction  and  multiple  knife-edge  diffraction  losses. 

A  detailed  description  of  the  SEKE  model  and  computer  listing 
of  the  program,  can  be  found  in  Reference  [1] . 

SEKE  was  designed  for  low  antenna  and  target  heights.  It 
has  been  tested  against  measurements  at  frequencies  ranging 
from  X-band  to  VHF  with  target  and  antenna  heights  between 
about  10  to  1000  m  [1].  When  used  for  target  and  antenna 
heights  outside  this  range  inaccuracies  and  failures  to  produce 
results  have  been  observed,  mainly  near  and  beyond  geometrical 
horizon.  These  inaccuracies  were  caused  by  the  spherical  earth 
diffraction  subroutine  (SPH35)  which  uses  a  series  due  to  Fock 
[2].  In  this  subroutine  the  values  of  the  Airy  function  in 
Fock' s  series  are  evaluated  by  a  fourth-order  polynominal  fit 
to  their  logarithms.  These  polynominals  are  tabulated  for 
normalized  effective  antenna  and  target  heights  of  up  to  100 
(i.e.  7000  m  at  VHF),  hence  limiting  the  applicability  of  SPH35 
to  low  effective  altitudes.  In  addition,  inaccuracies  in  the 
polynomial  fits  result  in  erroneous  values  for  the  propagation 
factor  at  some  combinations  of  antenna  and  target  heights  near 
the  optical  horizon,  when  many  terms  are  needed  in  Fock' s 
series . 

In  this  report,  we  present  a  more  accurate  method  for 
evaluating  the  Airy  function  in  the  complex  plane,  resulting  in 
the  extension  of  the  range  of  validity  of  SEKE.  The  Airy  func¬ 
tion  is  evaluated  using  the  power  series  expansion  near  the 
origin,  and  an  integral  representation  elsewhere.  The  computed 
values  of  the  Airy  function  agree  to  at  least  four  digits  with 
tabulated  values. 
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The  algorithm  that  evaluates  the  Airy  function  was  incor¬ 
porated  in  a  new  spherical  earth  diffraction  subroutine 
(SPH35N) .  It  was  found  that,  if  SEKE  uses  this  subroutine,  no 
problems  arise  for  normalized  heights  of  up  to  5000,  that  is, 
about  375  km  at  VHF  or  17  km  at  Ku  band. 

Some  background  information  is  given  in  Section  2  about 
Fock's  series,  and  about  the  Airy  integral  and  its  properties. 
Notation  is  also  introduced  in  this  section.  The  algorithm  for 
evaluating  the  Airy  function  is  described  in  Section  3  by  ex¬ 
plaining  the  power  series  expansion,  the  Gaussian  quadratures 
integration,  and  the  connection  formula.  Certain  computational 
difficulties  resulting  from  the  incorporation  of  this  Airy 
function  subroutine  in  the  spherical  earth  diffraction  algo¬ 
rithm  are  described,  and  the  method  used  to  overcome  them  is 
presented  in  Section  4.  Finally,  a  comparison  between  predic¬ 
tions  made  by  SEKE,  using  either  SPH35  or  SPH35N  are  provided 
in  Section  5.  A  flowchart  for  SPH35N  and  its  code  are  appended 
at  the  end  of  this  report. 

2 .  BACKGROUND 

2 . 1  SEKE 

Low-altitude  propagation  loss  is  influenced  by  atmospheric 
refraction  and  by  diffraction  and  multipath  (reflection)  from 
the  terrain  over  which  the  waves  travel.  As  a  result  SEKE  has 
to  make  two  main  decisions.  It  has  to  first  decide  whether  to 
use  a  line-of-sight  model  (multipath),  or  a  diffraction  model. 
If  the  program  decides  to  employ  a  diffraction  model,  a  second 
decision  must  be  made:  to  use  multiple  knife-edge  diffraction 
or  spherical  earth  diffraction. 

The  first  decision  is  based  on  the  clearance  of  the  direct 
ray  between  the  radar  and  the  target.  The  program  first  lo¬ 
cates  on  the  terrain  profile  the  highest  mask  of  the  minimum 
clearance  point,  M.  Let  A  be  the  clearance  of  the  direct  ray 
at  the  point  M  and  d-2  and  d2  be  the  distances  from  the  radar  to 
M  and  from  M  to  the  target  respectively.  The  Fresnel  clearance 
at  the  point  M  is  given  by  the  formula: 
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where  X  is  the  wavelength.  If  A/A0  >  1  SEKE  uses  multipath 
alone.  If  A/A0  <  1/2  SEKE  uses  the  diffraction  subroutines. 

For  the  intermediate  zone  where  1>  A/A0  >  1/2,  a  weighted 
average  of  the  multipath  loss  and  the  diffraction  loss  F^,  is 
used  (see  Reference  [1]  for  more  detail) .  When  A/A0  <  1  a 
second  decision  is  made  between  multiple  knife-edge  and  spheri¬ 
cal  earth  diffraction.  The  key  parameter  in  this  solution  is 
the  ratio  hj^/A0,  where  h^  is  the  height  of  the  highest  mask  or 
minimum  clearance  point  M,  measured  from  the  best-fit  line  on 
the  terrain  profile.  When  the  height  of  the  discrete  obstacles 
over  smooth  earth  are  large  relative  to  A0,  knife-edge  dif¬ 
fraction  should  dominate  over  spherical  earth  diffraction.  In 
SEKE  the  propagation  loss  is  approximated  by  knife-edge  dif¬ 
fraction  alone  when  hj^/A0  >  1/2,  by  spherical  earth  diffrac¬ 
tion  when  h^/A0<  1/4  and  a  weighted  average  of  the  two  is  used 
in  the  intermediate  region. 

2.2  Inaccuracies  in  SEKE 

Figure  1  shows  some  of  the  inaccuracies  in  SEKE  caused  by 
the  SPH35  routine.  The  plot  shows  the  two-way  propagation  fac¬ 
tor  over  a  smooth  conducting  spherical  earth  with  an  antenna 
height  of  10  m  and  a  target  height  of  8000  m  at  167  MHz.  At 
ranges  below  185  km,  SEKE  uses  the  geometric  optics  routine 
GEOSE .  After  that  point,  SEKE  attempts  to  use  SPH35,  but  that 
routine  fails.  SEKE  falls  back  in  GEOSE,  which  remains  fairly 
accurate  until  about  350  km.  Then  SEKE  starts  producing  incor¬ 
rect  results.  At  the  optical  horizon,  which  is  at  385  km,  SEKE 
is  still  attempting  to  use  SPH35,  which  still  fails,  but 
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Fig.  1.  Prediction  of  the  two  way  propagation  loss  as  a  function 
of  range  for  antenna  height  =  10  m  and  target  height  =  8000  m  at 
VHF  (167)  MHz).  The  solid  line  shows  the  prediction  according  to 
SPH35  and  the  dashed  line  according  to  SPH35N. 
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SEKE  now  uses  the  knife-edge  routine  LAPKE.  This  results  in 
erroneous  values.  At  a  range  of  620  km,  SPH35  returns  a  value 
acceptable  to  SEKE,  which  is,  however,  incorrect.  (The  dashed 
line  of  Figure  1  shows  the  values  returned  by  SPH35N) .  Further 
study  shows  that  the  problems  in  SPH35  are  associated  with  the 
computation  of  Airy  functions  in  that  routine. 

2.3  Spherical  Earth  Diffraction  Subroutine  and  Fock' s 


Series 


Fock  has  shown  [3],  that  the  propagation  factor,  F,  is  a 
function  of  normalized  antenna  height,  y,  normalized  target 
height,  z,  and  normalized  range,  x,  given  by  the  following  sum: 


(1) 


n=l 


F  is  the  propagation  factor,  defined  here  to  be  the  ratio  of 
the  electric  field  at  a  point  to  the  free  space  electric  field 
at  that  point. 


Man  +  exp(f  )v) 
Jn{U)  exp  (*)Ai'(an) 


(2) 


Ai (w)  is  the  Airy  integral: 


and  an  are  the  zeroes  of  the  Airy  integral. 


An  antenna  height  (ha)  and  a  target  height  (h^) ,  cor¬ 
respond  to  normalized  heights  of 


ha  ht 


(3) 
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respectively,  where  the  normalization  constant  is  given  by: 


(4) 


where  a  is  the  effective  radius  of  earth,  and  X  is  the 
wavelength.  The  range  normalization  factor  rQ  and  the  normal¬ 
ized  range  x  are  given  by: 


(5) 


SPH35,  the  algorithm  that  has  been  employed  to  evaluate 
the  propagation  factor  in  the  spherical  earth  diffraction 
region,  uses  the  first  n  terms  (n  >  35)  in  Fock's  series.  The 
convergence  criterion  is  that  the  contribution  of  each  one  of 
three  successive  terms  of  the  series  is  less  than 
8  =  0.0005.  If  this  criterion  is  not  met  while  considering  the 
first  35  terms,  SPH35  returns  a  message  that  it  "diverged". 

When  SPH35  "diverges"  SEKE  applies  the  geometrical  optics  sub¬ 
routine  if  the  point  is  visible,  otherwise  it  uses  the  knife- 
edge  diffraction  subroutine. 

SPH35  expresses  fn(u)  as  fn  (u)  =  exP  (^-n(u)  +  iM-n  <u)  )  > 
where  Xn  (u)  and  |ln(u)  are  fitted  by  fourth  order  polynomials  in 
u.  The  coefficients  of  the  polynomials  for  each  term  of  the 
series  are  tabulated  in  SPH35.  In  order  to  make  the  polynomial 
fit  more  accurate,  the  range  of  normalized  heights  (u)  is 
divided  in  four  separate  regions.  (See  Table  I).  fn(u)  =  u 
has  been  found  to  be  a  good  approximation  in  the  interval 
0  <  u  <  0.05.  For  the  regions,  0.05  <  u  <  0.2,  0.2  <  u  <  3, 

3  <  u  <  10,  and  10  <  u  <  100  four  different  polynomial  fits  are 
used.  No  attempt  is  made  in  this  original  version  of  SPH35  to 
obtain  accurate  polynomial  fits  to  A.n  and  |in  for  u  >  100.  In 
this  high  altitude  region,  SPH35  returns  coefficients  for  Xn 
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TABLE  I 


Approximations 

Used  in  SPH35  to  Evaluate  fn (u)  for  Different 
Ranges  of  Normalized 

Target  Height  u 

RANGE 

CALCULATION  of  fn (u) 

0  <  u  <  0.05 

fn  (u)  =  u 

0.05  <  u  <  0.2 

1st  fourth-order  polynomial  fit 

0.2  <  u  <  3 

2nd  fourth-order  polynomial  fit 

3  <  u  <  10 

3rd  fourth-order  polynomial  fit 

10  <  u  <  100 

4th  fourth-order  polynomial  fit 

u  >  100 

value  of  u  out  of  range  (however  SPH35 
returns  a  value  for  F  according  to  the 
4th  fourth-order  polynomial  fit) 
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and  ^in  according  to  the  fourth  polynomial  fit.  As  a  result, 
the  accuracy  of  SPH35  is  limited  to  normalized  heights  of  up  to 
100  (i.e.  7000  m  at  VHF,  or  317  m  at  Ku  band) . 

There  are  two  potential  problems  with  SPH35.  First,  there 
is  no  polynomial  fit  for  u  >  100.  Second,  even  for  lower 
altitudes,  inaccuracies  in  the  polynomial  fits  result  in 
erroneous  values  for  the  propagation  factor  or  failure  to  con¬ 
verge  for  some  combinations  of  antenna  and  target  heights  near 
the  optical  horizon.  These  discrepancies  occur  in  regions 
where  the  sum  of  the  series  is  much  smaller  than  the  largest 
term,  or  where  many  terms  are  required  for  convergence. 

In  order  to  test  the  accuracy  of  the  polynomial  fits,  the 
values  returned  by  SPH35  were  compared  with  the  values  returned 
by  another  subroutine  PROPSES.  PROPSES  calculates  the  Airy 
functions  by  numerically  integrating  the  differential  equations 
which  define  them.  This  technique  results  in  high  accuracy  but 
requires  long  computational  time. 

Table  II  shows  the  percent  error  in  SPH35.  Sample  runs 
are  presented  for  each  of  the  polynomial  fits.  It  is  seen  that 
the  fourth  polynomial  fit  is  inaccurate  even  far  into  the  dif¬ 
fraction  region.  The  remaining  fits  seem  to  work  adequately 
far  into  the  diffraction  region,  but  SPH35  sometimes  fails  to 
return  a  value  or  return  an  inaccurate  value  when  called  well 
inside  the  horizon. 

It  is  apparent  from  the  above  discussion  and  the  sample 
run  presented  that  the  accuracy  of  the  fourth  order  polynomial 
fits  for  fn(u)  is  not  always  sufficient  to  obtain  a  reasonable 
prediction  for  the  propagation  factor.  As  a  result,  there  is  a 
need  for  a  new  spherical  earth  diffraction  subroutine  where  the 
Airy  function  is  evaluated  using  a  more  accurate  method.  The 
remainder  of  this  report  describes  such  a  new  spherical  earth 
diffraction  subroutine. 
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TABLE  II 


SAMPLE 
OF  fn(u) 


RUNS  FOR  FINDING  THE  PERCENT  ERROR 
FOR  EQUAL  ANTENNA  AND  TARGET  HEIGHTS 


Method  for  ha  =  h^ 

calculating  (m) 

fn(u)  in 
SPH35 


%  Error  %  Error  at  a  range  of 

at  first  point  1.5  times  the 
where  spherical  optical  horizon 
earth 

diffraction 
is  used  by  SEKE 
(well  within  the 
optical  horizon) 


fn(u)  =  u 

3 

18.07 

-  0 

1st  fit 

10 

SPH35  failed 

0.15 

2nd  fit 

70 

-  0 

0.03 

3rd  fit 

600 

16.68 

0.72 

4th  fit 

5,  000 

SPH35  failed 

28.65 
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2.4  Properties  of  the  Airy  Function 

Algorithms  for  the  evaluation  of  the  Airy  function  with 
real  arguments  are  readily  available  [8],  However,  their  modi¬ 
fication  to  handle  complex  arguments  is  not  straightforward. 

The  power  series  expansion  can  be  used  only  in  a  small  region 
near  the  origin,  and  the  asymptotic  formulas  do  not  give  suffi¬ 
cient  accuracy  for  moderate  values  of  z  =  x  +  iy  [3] . 


The  Airy  function  Ai (z) 


satisfies  the  differential  equation 


The  two  sets  of  linear  independent  solutions  to  the  above  equa¬ 
tion  are  Ai(z),  Bi(z)  and  Ai[(z  exp  +  (2rci/3)  ]  where 


(6) 


The  functions  are  entire  and  so  have  convergent  power  series. 

Schulten  et  al  [4]  give  an  integral  representation  for  the 
*  Airy  function  whose  evaluation  by  a  Gaussian  quadrature  method 
requires  only  a  few  terms.  The  integral  representation  for 
Ai(z)  is  derived  from  an  expression  for  the  modified  Bessel 
function  of  the  second  kind  (z) .  (formula  6.627  of  Ref. [5] . 


If  we  set  b  =  1/3,  £  =  2  zJ/^  and  substitute: 

3 


(7) 


(8) 
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(9) 


equation  (7)  can  be  solved  for  Ai(z). 


Ai(z)  =  l-K -l/2z-l/4e-2/3z 


\argz\  <  ^  \z\  >  0 


p (x)  is  a  non-negative  exponentially  decreasing  function: 


(10) 


This  integral  representation  is  valid  in  the  sector 
I arg  z|  <  2n/3.  However,  there  exists  a  connection  formula 
that  transforms  a  point  outside  this  sector  to  a  weighted  sum 
of  two  linearly  independent  points  inside  it: 


Ai(z)  =  eTi/3Ai(ze~2Ti/3)  +  e~Tl/3Ai(ze2r,/3) 


(11) 


(See  formula  10.4.7  of  Ref.  [3])  . 

Even  though  p(x)  contains  the  Airy  function,  the  weights 
and  abcissas  for  the  Gaussian  quadrature  can  be  computed 
without  an  accurate  computation  of  the  Airy  function,  because 
the  moments  of  p(x)  can  be  evaluated  in  closed  form. 

3.  ALGORITHM  FOR  EVALUATING  THE  COMPLEX 
AIRY  FUNCTION 

As  mentioned  above  there  exists  a  power  series  expansion 
for  the  Airy  function.  This  series  method  is  used  to  evaluate 
the  value  of  the  Airy  function  close  to  the  origin.  For  large 
values  of  z  a  Gaussian  quadratures  method  is  implemented  to 
evaluate  the  integral  representation  above.  These  two  numeri¬ 
cal  methods  are  used  in  the  part  of  the  complex  plane  where 
I arg  z|  <  2n/3.  For  the  remaining  part  of  the  plane  the  con¬ 
nection  formula  is  used. 

It  was  found  by  looking  at  the  values  returned  by  these 
two  different  methods  (power  series  and  Gaussian  quadratures) 
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lm(z) 


INTEGRAL  REPRESENTATION 
POWER  SERIES 
CONNECTION 


Re(z) 


Fig.  2.  Regions  on  the  complex  plane  where  the  power  series 
expansion,  the  Gaussian  quadratures  method  to  solve  the 
integral  representation  and  the  connection  formula  are  each 
used. 
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on  certain  radials  of  the  right  half-plan  (x>0)  that  the  two 
methods  returned  the  same  value  for  Ai(z),  within  four  sig¬ 
nificant  digits,  when  |z|  =  2.  As  a  result,  the  algorithm  was 
implemented  so  that  if  z  lies  in  the  right  half-plane  and 
|z|  £  2,  then  Ai  (z)  is  evaluated  by  the  power  series  expan¬ 
sion.  If  |z|  >2  (and  x>  0)  then  Ai(z)  is  evaluated  by  the 
Gaussian  quadratures  method.  When  x  <  0  (left-hand  plane) ,  it 
was  found  that  the  two  methods  overlapped  at  | z |  «  4 .  As  a 
result  if  |z|  >4  and  x  <  0,  then  the  power  series  expansion 
algorithm  is  used.  If  |z|  >4  and  x  <  0,  the  Gaussian 
quadratures  method  is  applied.  Figure  2  defines  the  region  on 
the  complex  plane  where  the  power  series  expansion,  the  Gauss¬ 
ian  quadratures  method  and  the  connection  formula  are  used. 

3.1  Power  (Taylor)  Series  Expansion 

The  Airy  function  is  entire,  so  a  convergent  Taylor  series 
representation  exists  for  it  (See  Eq.  10.4.2  of  Ref.  [3]). 

This  series  converges  very  fast  for  small  values  of  z.  When  a 
large  value  of  z  is  used,  the  series  converges  very  slowly,  and 
inaccuracies  occur  due  to  large  cancellations.  The  Taylor 
series  expansion  is  given  by: 


Ai(z)  =  a  ■  h(z )  -  0  ■  g(z ) 


(12) 


where : 


a  =  .Az(O)  — 


Bi(0) 

v/3 


3-2/3 

W/3) 


0.3550280538 


(13) 
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(14) 


(15) 


(16) 


3‘(7  +  i)  =1 
3  o 

3‘(7  +  1)  =  (37  +  1)(37  +  4)  •  •  •  (3T  +  3k  - 2) 

o  k 


(17) 


(7  arbitrary,  k=l ,2,  •  •  •) 


Equation  (12)  is  implemented  in  the  algorithm  that 
evaluates  the  Airy  Function  in  the  power  series  region  (see 
Figure  2 ) . 

In  order  to  reduce  the  growth  of  the  round-off  error,  the 
series  are  evaluated  term  by  term: 


m 


Ai(z)  =  £  an(z) 


n=  1 

for  an  m  such  that  am(z)  <  10-10,  where 
an(z)  =  a  hn(z)  -  p  gn(z) 

and  hn(z)  and  gn(z)  are  the  nth  terms  in  the  expressions  for 


g (z)  and  h (z) . 

3.2  N-Point  Gaussian  Quadratures  Approximation 

As  mentioned  above  Ai(z)  is  evaluated  for  large  z  by  the 
generalized  Gaussian  quadratures  method.  The  challenging  part 
with  this  method  is  to  find  a  function  p(x)  so  that  the  in- 
tegrable  singularities  are  removed  from  the  Airy  integral. 

Given  p(x)  and  given  the  number  of  terms  N  that  should  be  used, 
one  can  find  a  set  of  weights  w-^  and  abscissas  x-[  such  that  the 
approximation : 
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(18) 


[a  - 
/  p(x)p(x)dx  ~  22w;p(xi) 

Jb  «=i 


is  exact,  if  p(x)  is  a  polynomial  of  degree  less  than  2N. 
Schulten  et  al.[4]  give  an  integral  representation  for  Ai(z), 
whose  evaluation  by  a  Gaussian  quadrature  method  requires  only  a 
few  terms : 


Ai(z)  =  lr-i/2*-i/4e-C  r  fW-dx 
2  Jo  x  +  C 


\z\  >  0  and  |  arg£|  < 


where 


(19) 


p(x)  =  7T -'t22-u'63-2'3x-2'3e-*Ai 


'3x^2/3' 


and  £  =  -z3/2 


(20) 


The  integral  portion  of  equation  (19)  can  be  approximated 
by  the  quadrature  formula: 


Jr  °o 
' 

o 


p(j) 

C  +  X 


n 


w« 
C  + 


(21) 


The  weights  w-^  and  the  abscissas  were  found  by  im¬ 
plementing  the  procedure  described  by  Press  [6]. 

The  "scalar  product  of  two  functions  f  and  g  over  a  weight 
function  p(x)",  is  defined  as: 


< 


rb 

f\ 9  >=  /  p(x)f(x)g(x)dx 

J  a 


(22) 


A  set  of  orthogonal  polynomials  that  includes  exactly  one 
polynomial  of  order  j,  called  Pj (x)  j=0,  1,  2,  ...  is 

needed  to  find  the  weights  and  abscissas. 
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This  set  of  polynomials  can  be  constructed  by  the  following 
recurrence  relation: 


Po(x)  =  1 


(23) 


(the  second  term  is  omitted  when  i=0) ,  and 


(24) 


The  integral  of  equation  (24)  can  be  readily  evaluated  by 
observing  that  p(x)  is  a  solution  to  the  Stieltjes  moment  prob¬ 
lem,  whose  moments  (i^  can  be  explicitly  evaluated:  (Formula 

6.621.3,  [5]) 


k  =  0,1,2,  ... 


(25) 


as  a  result  <xp^|p^>  becomes  a  sum  of  M-k's.  Once  the  abscissas 


X]_,  X2  ....  xjj  (i.e.  the  N  zeroes  of  p^j(x))  are  known,  the 
weights  w-^  can  then  be  found. 


Press  [6]  presents  a  simple  method  to  find  the  wj_'s.  A 
new  sequence  of  polynomials  cp(x)  is  constructed,  by  the  fol¬ 
lowing  recurrence: 


<Po(x)  =  0 

=  p 
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where  p'  is  the  derivative  of  p  (x) . 
*1  1 


The  weights  of  an  N-point  Gaussian  quadrature  are  given  by 
the  relation: 


il\  = 


V>n(x  >) 

PN(xi) 


i=l,2,  •••,  N 


(27) 


The  procedure  described  above  was  implemented  to  find  the 
weights  w-^  and  abscissas  x-^  for  the  N-point  Gaussian  quadrature 
approximation.  When  the  Airy  function  is  evaluated  by  this 
method  it  is  dependent  on  N,  the  number  of  terms  used  in  the 
summation: 


Ai(z,  N) 


1  ^  N 
1  £ 


1=0 


C  +  xi 


where  £  =  - 2 2//3 


(28) 


We  calculated  Ai  (z,  N)  (for  |z|  >4)  for  N  =1,  2  ...20.  It 
was  found  that  the  values  of  Ai(z,  N)  were  almost  the  same  for 
values  of  N  near  10.  For  larger  or  smaller  N,  the  values  of  Ai 
(z,  N)  were  quite  different  indicating  truncation  or  round-off 
error.  It  was  thus  decided,  that  N=10,  is  a  good  choice  for 
the  number  of  terms  that  should  be  used  in  the  Gaussian  quadra¬ 
tures  approximation.  Table  III  gives  a  list  of  the  weights  w-j_ 
and  abscissas  x^  for  the  10-term  Gaussian  Quadrature  integra¬ 
tion  for  the  Airy  function. 

4 .  ALGORITHM  FOR  EVALUATING  THE  PROPAGATION  FACTOR  IN  THE 
SPHERICAL  EARTH  DIFFRACTION  REGION  (SPH35N) 

4.1  Computational  Difficulties  and  Solutions 

The  above  algorithm  for  evaluating  the  Airy  function  was 
implemented  in  a  new  spherical  earth  diffraction  subroutine 
(SPH35N) .  However,  the  incorporation  of  this  algorithm  to 
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TABLE  III 


10 -TERM  GAUSSIAN  QUADRATURES  INTEGRATION  FOR  THE  AIRY  FUNCTION 


Abscissas  Weights 


i  xi  wj_ 


1 

1.408308107197377E+01 

2. 677084371247434E-14 

2 

1. 02148854806031 5E+01 

6. 63 67 68 68 8 1 75 8 7 0E- 11 

3 

7 .441601846833691E+00 

1.7584 05638 61 9854E-0 8 

4 

5.307094307915284E+00 

1 .37 123 914 8 97 684 8E- 06 

5 

3. 634013504378772 E+ 00 

4 .435096659959217E-05 

6 

2. 331065231384 954E+00 

7 .155501 0754 31 907E-04 

7 

1. 34479708313994 5E+00 

6. 488956601264211 E- 03 

8 

6. 418885840366331 E- 01 

3.644041585109798E-02 

9 

2 .01 0034 600 90571 8E-01 

1 .439979241604145E-01 

10 

8.059435921534400E-03 

8. 1231 14 13423598 0E-01 
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compute  fn(u)  in  Fock's  series,  is  not  straightforward.  A 
problem  arises  in  the  evaluation  of  Ai (w) ,  where 
w  =  an  +  exp(7ti/3)  u,  by  the  Gaussian  quadrature  method.  As 
I w |  becomes  large,  in  the  region  |arg  wj  >  7t/ 3 ,  the  expression 
exp  (-2/3  w3/2)#  could  overflow.  However,  looking  at  Fock's 
series  the  term  exp (1/2 (V3  +  i)  anx)  can  be  used  to  partially 
cancel  a  large  value  for  exp (-2/3  w^/2) .  The  subroutine  that 
evaluates  the  Airy  function  returns 

aT(w)=  Ai (w)  exp  (2/3  w^/2) f  instead  of  Ai (w) .  The  Fock's 
series  equation  implemented  in  SPH35N,  then  becomes: 


OO  _  _  r  /  J  \  ■ 

F(x ,  y,  z)  =  2^/t tx  h(y)Mz) exp  ~1>  ~  C  +  (  ^(V3  +  i)anxj 

n= 1 


(29) 


where 


—  Ai  (an  +  o  o 

"  -kmfcj 1  •  *  =  f3/;  c = 1^> 


The  connection  formula  used  when  Ai(z)  is  evaluated  instead  of 
Ai  (z) ,  becomes : 


Ai(z)  =  e”l3Tt  (ze-ln>3)  +  f^Ai  (zS’i/3) 

where  Ai(w)  =  Ai(w)c 2^3u 

SPH35N  considers  as  many  terms  in  Eq.  (29)  as  are  needed 
to  get  two  consecutive  terms  that  contribute,  in  absolute 
value,  less  than  0.0005.  If  this  does  not  occur  within  the 
first  35  terms,  SPH35N  is  considered  to  have  "diverged".  One 
more  check  is  performed  in  SPH35N  to  make  sure  that  the  com¬ 
putation  is  accurate.  If  any  one  of  the  terms  in  equation  (29) 
contributes  more  than  10,000,  SPH35N  is  interrupted.  A  maximum 
contribution  of  10,000  per  term  was  picked  since  Ai (w)  is 
evaluated  with  an  accuracy  of  at  least  four  significant  digits. 
This  criterion  ensures  that  SPH35N  does  not  produce  wholly 
spurious  results  because  of  cancellation. 
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Empirically,  10,000  seemed  to  work  very  well.  SEKE  does 
try,  for  some  few  combinations  of  antenna  and  target  heights, 
to  use  SPH35N  in  the  multipath  region  far  above  the  optical 
horizon.  When  this  happens,  a  term  in  SPH35N  contributes  more 
than  10,000,  making  it  impossible,  given  the  accuracy  of  the 
Airy  function,  to  cancel  its  contribution.  Had  this  check  not 
been  included,  SPH35N  would  have  incorrectly  concluded  that  the 
series  converged.  With  this  check,  SEKE  uses  geometric  optics, 
which  is  accurate  in  this  region. 

In  certain  other  cases,  where  a  large  obstacle  masks  the 
target,  SEKE  tries  to  use  a  combination  of  spherical  earth  dif¬ 
fraction  and  knife-edge  diffraction.  In  some  of  these  cases, 
even  though  the  target  is  masked,  it  can  be  well  above 
the  horizon  for  the  smooth  fit  to  the  terrain  used  by  SEKE  to 
determine  the  effective  heights  for  the  spherical  Earth 
routine;  this  sometime  results  in  SPH35N  "diverging".  When 
this  happens,  SEKE  is  forced  to  use  only  knife-edge  diffraction 
as  expected,  since  the  dominant  effect  is  the  obstacle. 

4.2  Range  of  Validity  of  SPH35N  when  Called  from  SEKE 

The  range  of  validity  of  SPH35  in  target-to-antenna  dis¬ 
tance  x,  was  investigated  for  a  fixed  normalized  target  height, 

y,  and  antenna  height,  z.  It  was  found  that  the  range  where 
SPH35N  broke  down  was  always  in  a  region  where  SEKE  would  use 
geometrical  optics,  since  the  minimum  clearance  was  found  to  be 
greater  than  a  Fresnel  clearance  for  any  combination  of  y  and 

z,  where  y  <  5000  and  z  <  5000  (i.e.  heights  of  less  than  375 
km  at  VHF  or  17  km  at  Ku  band) .  As  mentioned  above,  in  prac¬ 
tice  SEKE  occasionally  calls  the  sphere  diffraction  program  in 
a  region  where  it  fails  to  return  an  answer.  SEKE  then  either 
uses  the  geometric  optics  calculation  when  the  point  is  in  the 
multipath  region  or  the  knife-edge  diffraction  calculations 
when  the  point  is  masked.  It  has  been  found  empirically  that 
SEKE  makes  an  accurate  decision  in  such  cases. 


20 


5.  EVALUATION  OF  SPH35N 

SPH35N  was  run  for  certain  cases  where  SPH35  would  not 
work  correctly  either  because  y  or  z  were  greater  than  100,  or 
because  of  cancellation  error  when  y  and  z  are  less  than  100. 
Two  examples  are  shown  in  Figures  3  and  4 .  Figure  3  presents 
the  two-way  propagation  factor  F^  in  dB  with  respect  to  range 
for  antenna  heights  ha  =  10  m  (z  =  0.14)  and  target 
height  h^  =  2.000  m  (y  =  28.57)  at  VHF.  The  solid  line  gives 
the  prediction  using  SPH35  and  the  dashed  line  the  prediction 
using  SPH35N.  The  spike  in  the  prediction  disappears  when 
SPH35N  is  used.  Figure  4  presents  another  plot  of  F^  (dB)  with 
respect  to  range  for  antenna  height  ha  =  6000  m  and  target 
height  ht  =  10  m.  In  this  case,  the  plot  resulting  by  the  use 
of  SPH35  shows  most  of  the  problems  created  by  the  inaccurate 
evaluation  of  the  Airy  function  in  the  spherical  earth  diffrac¬ 
tion  subroutine.  At  around  275  km,  SEKE  tries  to  use  spherical 
earth  diffraction,  but  SPH35  "diverges",  so  SEKE  uses  the  geom¬ 
etrical  optics  subroutine  up  to  the  optical  horizon.  Beyond 
the  optical  horizon  knife-edge  diffraction  is  used,  returning  a 
constant  value  of  F^  =  -6  dB.  At  around  350  km  SPH35  "con¬ 
verges"  and  returns  a  value.  When  SPH35N  is  employed,  the 
resulting  curve  is  smooth.  Problems  as  the  ones  encountered 
when  using  SPH35  do  not  occur,  since  SPH35N  is  able  to  return  a 
value  when  called  by  SEKE,  as  expected  from  the  specified  range 
of  validity  of  SPH35N . 

6.  CONCLUSIONS 

We  have  presented  a  subroutine  that  evaluates  the  propaga¬ 
tion  factor  in  the  spherical  earth  diffraction  region  by  apply¬ 
ing  Fock's  series  using  an  accurate  and  efficient  algorithm  to 
evaluate  the  complex  Airy  function. 

The  Airy  function  is  evaluated  by  the  power  series  expan¬ 
sion  for  |z|  close  to  the  origin,  by  the  10  term  Gaussian 
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RANGE  (km) 


Fig.  3.  Prediction  of  the  two  way  propagation  loss  as  a 
function  of  range  for  antenna  height  =  10  m  and  target  height 
=  2000  m  at  VHF  (167  MHz) .  The  solid  line  shows  the  prediction 
according  to  SPH35  and  the  dashed  line  according  to  SPH35N. 
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Fig.  4.  Prediction  of  the  two  way  propagation  loss  as  a 
function  of  range  for  antenna  height  —  10  m  and  target  height 
=  6000  m  at  VHF  (167  MHz) .  The  solid  line  shows  the  prediction 
according  to  SPH35  and  the  dashed  line  according  to  SPH35N . 
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quadrature  approximation  for  large  |z|  and  the  connection  for¬ 
mula  for  |  arg  z|  <  (2/3)  7t  i.e.  in  the  region  where  the  integral 
representation  is  invalid.  The  implemented  Airy  function 
checks  with  Airy  tables  within  four  significant  digits.  Incor¬ 
porating  this  subroutine  in  SEKE  eliminated  inaccuracies  caused 
by  the  old  spherical  earth  diffraction  subroutine  (SPH35) ,  that 
uses  a  fourth-order  polynomial  fit  to  approximate  the  Airy 
function.  The  new  subroutine,  SPH35N,  was  found  to  be  equally 
fast  as  SPH35  and  adequately  accurate  for  normalized  heights  of 
less  than  5000,  so  that  it  performs  as  desired  when  called  from 
SEKE. 
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Listing  of  SPH35N 


C***************************************************************** 


C  * 
C  SPHERICAL  EARTH  DIFFRACTION  LOSS  SUBROUTINE  * 
C  APPLYING  THE  ANALYSIS  DESCRIBED  IN  * 
C  FOCK ,  V .  A  . , ELECTROMAGNETIC  DIFFRACTION  AND  PROPAGATION  PROBLEMS , * 
C  OXFORD :PERGAMMON  PRESS,  LTD  ,  1965  * 
C  * 
C  GEORGE  H.  POLYCHRONOPOULOS  1  /87  * 
C  * 


c***************************************************************** 


SUBROUTINE  SPH35N  (RNG ,Z1 ,Z2, WAVE, AE.F) 

REAL  RNG, Z1,Z2, WAVE, AE.F 

REAL  *8  X , Y ,Z , COEF1 ,COEF, SENSI , ATER, APTERM ,DO , HO 
COMPLEX* 16  C0EF2 , EPI , FNY , FNZ , TERM , ETERM , SUM , FAIRY , AY , AZ , 
+  ZETAY , ZETAZ 

INTEGER  N .FLAG 
DIMENSION  A(35) ,DA(35) 


DATA  PI  /3. 141592/ 

C  NEGATIVE  OF  THE  ZEROES  OF  THE  AIRY  FUNCTION 


DATA  A  /-2 

.33810 

,  -4.08794, 

-5 

.52055, 

-6 

.78670 

,  -7.94413, 

+ 

-9 

.02265 

,-10.04017, 

-11 

.00852, 

-11 

.93601 

,-12.82877, 

+ 

-13 

.69148 

,-14.52782, 

-15 

.34075, 

-16 

. 13268 

,-16.90563, 

+ 

-17 

.66130 

,-18.40113, 

-19 

.12638, 

-19 

.83812 

,-20.53733, 

+ 

-21 

.22482 

,-21.90136, 

-22 

.56761, 

-23 

.22416 

,-23.87156, 

+ 

-24 

.51030 

,-25.14082, 

-25 

.76353, 

-26 

.37880 

,-26.98698, 

+ 

-27 

.58838 

,-28.18330, 

-28 

.77200, 

-29 

.35475 

,-29.93176/ 

THE 

DERIVATIVE  OF  THE  AIRY  FUNCTION  EVALUATED  AT  THE  ZEROES 

DATA  DA/O . 

70121, 

-0.80311, 

0. 

86520, 

-0. 

91085, 

0.94733, 

+ 

-0. 

97792, 

1.00437, 

-1. 

02773, 

1. 

04872, 

-1.06779, 

+ 

1. 

08530, 

-1.10150, 

1. 

11659, 

-1. 

13073, 

1.14403, 

+ 

-1. 

15660, 

1.16853, 

-1. 

17988, 

1. 

19070, 

-1.20106, 

+ 

1. 

21098, 

-1.22052, 

1. 

22970, 

-1. 

23854, 

1.24708, 

+ 

-1. 

25534, 

1.26334, 

-1. 

27109, 

1. 

27861, 

-1.28592, 

+ 

1. 

29302, 

-1.29994, 

1. 

30667, 

-1. 

31324, 

1.31965/ 

C  SET  UNDERFLOW  CONDITION  TO  NO  PRINTOUT 
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CALL  ERRSET(208, 0,-1, 1,1,0) 


C  FIND  THE  NORMALIZED  EQUIVALENT  OF  RNG,Z1,Z2 
D0=( AE*AE*WAVE/ 1E5/PI) ** . 33333333 
H0= ( AE*WAVE*WAVE/ 10/PI/PI) ** . 33333333/2 
X=RNG/DO 
Y=Z2/H0 
Z=Z1/H0 

C  EXP  (PI*I/3) 

EPI  =  (  0 . 5000001 8 12D+00 ,  0 . 8660252991D+00) 

C  C0EF1  =  2*SqRT(PI) 

C0EF1  =  0 . 3544908524D+01 
C  C0EF2  =  1/2*(SQRT(3)+I) 

C0EF2  =  (  0 . 8660254478D+00 ,  0 . 5000000000D+00) 

COEF  =  C0EF1  *  DSQRT  (X) 

SENSI=  0.0005 

C  INITIALIZE 

FLAG  =  0 
N  =  1 

APTERM  =  10000. 

SUM  =  (0. ,0.) 

C  COMPUTE  35  TERMS  IN  FOLK'S  SERIES  OR  UNTIL  THE  CONTRIBUTION  OF 
C  OF  TWO  SUCCESSIVE  TERMS  IS  LESS  THAN  0.0001  FOR  EACH  ONE 

2  IF  ((N.EQ.36) .OR. (FLAG.EQ.-2) .OR. (FLAG.EQ.-l))  GOTO  1 

AY  =  A(N)+EPI*Y 

FNY  =  FAIRY (AY )/(EPI*DA(N) ) 

AZ  =  A(N)+EPI*Z 

FNZ  =  FAIRY ( AZ) / (EPI*DA (N ) ) 

C  THE  VALUES  RETURNED  BY  FAIRY  ARE  AI(W)*EXP(ZETAW) ,  WHERE 
C  ZETAW  =  2./3.*(W**(3./2.)) 

ZETAY  =  2./3.*(AY*CDSQRT(AY)) 

ZETAZ  =  2./3.*(AZ*CDSQRT(AZ)) 

ETERM  =  -ZETAY-ZETAZ+ (C0EF2+A (N) *X) 

TERM  =  FNY*FNZ*CDEXP (ETERM) 

ATERM  =  CDABS  (TERM) 

IF  ( ATERM. LE. 10000)  GOTO  3 

C  FORCE  THE  PROGRAM  TO  TERMINATE  IF  ONE  TERM  CONTRIBUTES 

C  MORE  THAN  1000 

FLAG  =-2 

3  CONTINUE 

IF  ((APTERM. GE.SENSI) .OR. (ATERM. GE.SENSI))  GOTO  4 
FLAG  =-l 
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4  CONTINUE 
APTERM  =  ATERM 
SUM  =  SUM  +  TERM 
N  =  N+l 

GOTO  2 

1  CONTINUE 

F  =  COEF  *CDABS(SUM) 

C  F=-l  WHEN  35  TERM  ARE  NOT  ENOUGH  I.E  SPH35N  "DIVERGED" 
IF  (FLAG. NE. 0)  GOTO  5 
F  =  -1. 

5  CONTINUE 

C  F=-2  WHEN  ONE  TERM  CONTRIBUTES  MORE  THAN  10000 
IF  (FLAG . NE . -2)  GOTO  6 
F  =  -2 . 

6  CONTINUE 
RETURN 
END 


C************* ************************** ******** ********** ******** 


C  FUNCTION  TO  CALCULATE  AIRY  FUNCTIONS  * 

C  IN  THE  COMPLEX  PLANE  * 
C  THE  POWER  SERIES  AND  GAUSSIAN  QUADRATURE  METHODS  ARE  USED  IN  * 
C  THE  UPPER  HALF  PLANE.  THE  PROPERTY  OF  COMPLEX  CONJUGATE  OF  * 
C  AIRY  FUNCTIONS  IS  USED  TO  CALCULATE  THE  AIRY  FUNCTION  IN  THE  * 
C  LOWER  PLANE.  * 
C  * 
C  *  THIS  FUNCTION  RETURNS  AI (Z)*EXP(2* (Z**3/2)/2)  * 
C  * 


C************* ***************************** *********************** 

COMPLEX  FUNCTION  FAIRY*16  (Z) 


REAL* 8  RZ,IZ,TANZ 

C0MPLEX*16  Z, AIR, AI ,EPI , EPINT ,EPIN ,EPIT , AZETA 


C  CALL  ERRSET  TO  SUPPRESS  MESSAGES  WHEN  ONE  OF  THE  TERMS  UNDERFLOWS 
C  IN  THE  CONNECTION  FORMULA 

CALL  ERRSET  (208,0,-1,1,1,0) 


C  ASSIGN  THE  EXP  NEEDED  FOR  THE  CONNECTION  FORMULA 
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c 


c 

c 

c 


EXP  (PI*I/3) 

EPI  =  (  0 . 5000001812D+00 ,  0 . 8660252991D+00) 

EXP  (-2*PI*I/3) 

EPINT=  (-0.4999996375D+00,  -0 . 8660256131D+00) 
EXP  (-PI*I/3) 

EPIN  =  (  0 . 5000001812D+00 ,  -0 . 8660252991D+00) 

EXP  (-4*PI*I/3)  =  EXP  (2*PI*I/3) 

EPIT  =  (-0 . 4999996375D+00 ,  0 . 8660256131D+00) 

RZ  =  DREAL  (Z) 

IZ  =  DIMAG  (Z) 


C  IF  WE  ARE  IN  THE  CONNECTING  FORMULA  WEDGE  2PI/3<=ARG  Z<=4PI/3 
C  APPLY  MAPPING  FORMULA  TO  GET  OUT  OF  THE  WEDGE 

AZETA  =  CDEXP  (4 . *(Z*CDSQRT(Z))/3 . ) 

TANZ  =  IZ/RZ 

C  AI(Z)  RETURNS  AIRY(Z)*EXP(2/3*Z**2/3) 

IF  ((RZ.LT.O) .AND. (TANZ. LE. 1.7032) .AND. 

+  (TANZ. GE. -1.7032))  GOTO  1 

AIR  =  AI  (Z) 

GOTO  2 

C  CONNECTION  FORMULA  FOR  AI(Z) 

1  AIR  =  EPI*AI(Z*EPINT)*AZETA+EPIN*AI(Z*EPIT) 

2  CONTINUE 
FAIRY  =  AIR 
RETURN 

END 


COMPLEX  FUNCTION  AI*8  (AZ) 

C0MPLEX*16  AZ, POWER, GQA 
REAL*8  MAZ.RAZ 

C  SET  UNDERFLOW  CONDITION  TO  NO  PRINTOUT 

CALL  ERRSET(208, 0,-1, 1,1,0) 

MAZ  =  CDABS  (AZ) 

RAZ  =  DREAL  (AZ) 

C  IF  2PI/3  <=ARG  ZA  <=  PI/2  OR  4PI/3  <=  ARG  ZA  <=  3PI/2 
IF  (RAZ  .LE .0)  GOTO  3 
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IF  (MAZ.LE.2)  GOTO  7 
AI  =  GQA  (AZ) 

GOTO  8 

7  AI  =  POWER  (AZ) 

8  CONTINUE 
GOTO  4 

3  IF  (MAZ  .  LE.  4)  GOTO  5 

AI  =  GQA  (AZ) 

GOTO  6 

5  AI  =  POWER  (AZ) 

6  CONTINUE 

4  CONTINUE 

C  IF  IN  THE  RIGHT  HALF  PLANE  (  REAL(ZA)  >=  0) 
RETURN 
END 


C*** *********************************************** *************** 
C  POWER  SERIES  OF  AIRY  FUNCTION  + 
C  POWER  SERIES  ARE  COMPUTED  ITERATIVELY  UNTIL  THE  ADDITIONAL  TERM* 
C  CONTRIBUTES  LESS  THAN  1*E-10  * 
C  * 
C  THE  FORMULA  USED  WAS  TAKEN  FROM  J.C. MILLER  "THE  AIRY  INTEGRAL",* 
C  BRITISH  ASSOC.  ADV .  SCI.,  MATH  TABLES  VOLUME  B.1946  B17  * 

C* ******************************************************** ******** 


COMPLEX  FUNCTION  POWER* 16  (Z) 


REAL*8  ALPHA , BETA , DEN 1 , DEN2 , NUM1 , NUM2 , ATERM . FACTOR 
+  IN31 ,  IN32 

COMPLEX* 16  Z , AIRY , TERM1 , TERM2 , TERM , Y1 , Y2 , P 1 , P2 , EZETA 
INTEGER+4  P0WER1 , P0WER2 , FP 1 , FP2 , N 

C  SET  UNDERFLOW  CONDITION  TO  NO  PRINTOUT 

CALL  ERRSET(208, 0,-1, 1,1,0) 

ALPHA  =  0.355028053887817 
BETA  =  0.258819403792807 
ATERM  =  100000. 

EZETA  =  CDEXP  (2.*(Z*CDSQRT(Z))/3. ) 

N  =  1 

AIRY  =  ALPHA-(BETA+Z) 
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P0WER1=  3 
P0WER2=  4 

C  CONTINUE  ADDING  TERMS  UNTIL  THE  CONTIBUTION  OF  THE  LAST  TERM  IS 
C  LESS  THAN  IE-10 

2  IF  (ATERM.LE.1E-10)  GOTO  1 

Pl=Z**POWERl 
P2=Z**P0WER2 
FP1  =  POWER1 
FP2  =  P0WER2 
DEN1=  FACTOR(FPl) 

DEN2=  FACT0R(FP2) 

NUM1=  IN31 (N) 

NUM2=  IN32(N) 

TERM1  =  ALPHA* (NUM1 /DENI *P1 ) 

TERM2  =  BETA  * (NUM2/DEN2*P2) 

TERM  =  TERM1-  TERM2 
AIRY  =  AIRY  +  TERM 
POWERl=  P0WER1+3 
P0WER2=  POWER1+1 
N  =  N+l 

ATERM  =  CDABS(TERM) 

GOTO  2 

1  CONTINUE 

POWER  =  AIRY+EZETA 
RETURN 
END 
C 

o ***************************************************************** 

C - FUNCTION  THAT  RETURNS  THE  FACTORIAL  OF  AN  INTEGER - * 

C* ********************** ******************************************* 

c 

REAL  FUNCTION  FACT0R*8(N) 

C 

C  SET  UNDERFLOW  CONDITION  TO  NO  PRINTOUT 
CALL  ERRSET(208 , 0 , -1 , 1 , 1 p  0) 

IF  ((N.EQ.O)  .OR.  (N.EQ.l))  GOTO  10 
FACTOR  =  1. 

30  IF  (N.EQ.l)  GOTO  20 

FACTOR  =  FACTORS 
N  =  N-l 
GOTO  30 

20  RETURN 

10  FACT0R=1 . 

RETURN 
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END 


C 

C 

o ********* ******** ************** *************** ****************** 

C - FUNCTION  THAT  EVALUATES  1*4*7 _ 3*(N-l)  +  l - * 

C* ********************************************************** ****** 

c 

REAL  FUNCTION  1*31*8(1) 

C 

C  SET  UNDERFLOW  CONDITION  TO  NO  PRINTOUT 

CALL  ERRSET(208, 0,-1, 1,1,0) 

IN31  =1 

IF  (N.EQ.l)  GOTO  50 
DO  40  1=2, N 

IN31  =  IN31*(3*(I-1)+1) 

40  CONTINUE 

50  RETURN 

END 
C 
C 

c* ********* ******************************************************* 

C - FUNCTION  THAT  EVALUATES  2*5*8 _ 3*(N-l)+2 - * 

C* ******************************************************** ******** 

c 

REAL  FUNCTION  IN32*8(N) 

C 

C  SET  UNDERFLOW  CONDITION  TO  NO  PRINTOUT 

CALL  ERRSETC208, 0,-1, 1,1,0) 

IN32  =2 

IF  (N.EQ.l)  GOTO  60 
DO  70  1=2, N 

IN32  =  IN32*(3*(I-l)+2) 

70  CONTINUE 

60  RETURN 

END 


C* ******************************************************** ******** 


C  PROCEDURE  TO  CALCULATE  THE  AIRY  FUNCTION  OF  * 
C  A  COMPLEX  ARGUMENT  BY  THE  GAUSSIAN  QUADRATURES  METHOD  * 
C  AI(Z)  =  l/2*(PI**(-l/2))*(Z**(-l/4))*EXP(ZETA)*SUM  OVER  N  * 
C  W(I)/(1+X(I)/ZETA) )  * 
C  * 
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C  IMPLEMENTED  AS  IN  "AN  ALGORITHM  FOR  THE  EVALUATION  OF  * 
C  OF  COMPLEX  AIRY  FUNCTIONS" ,  JOURNAL  OF  COMPUTATIONAL  * 
C  PHYSICS  31,  60-75  (1979)  * 
C  * 
C  THIS  FUNCTION  RETURNS  AI (Z)*EXP (ZETA)  * 
C  * 
C  WEIGHTS  AND  ZEROES  WERE  CALCULATED  BY  GQAIRY . PLI  * 


C* ********************************************************** ****** 


COHPLEX  FUNCTION  GQA*16  (Z) 

REAL  *8  ZEROES  (10) .VEIGHT(IO) 
COHPLEX  *16  AIR, Z, SUM, ZETA 
INTEGER  N 

DATA  DSqRPI  /O. 564189584/ 

C  SET  UNDERFLOW  CONDITION  TO  NO  PRINTOUT 
CALL  ERRSET(208, 0,-1, 1,1,0) 


C  INITIALIZE  (ASSIGN  VALUES  TO  WEIGHTS  AND  X-INTERCEPT) 

ZER0ES(1)=  1.408308107197377E+01 
ZEROES (2) =  1 . 021488548060315E+01 
ZER0ES(3)=  7.441601846833691E+00 
ZEROES (4) =  5 . 307094307915284E+00 
ZER0ES(5)=  3 . 634013504378772E+00 
ZER0ES(6)=  2 . 331065231384954E+00 
ZEROES (7) =  1 . 344797083139945E+00 
ZEROES (8) =  6 . 418885840366331E-01 
ZEROES ( 9) =  2 . 0100346009057 18E-01 
ZEROES (10)=  8 . 05943592 1534400E-03 
WEIGHT(l)=  2 . 677084371247434E-14 
WEIGHT(2)=  6 . 636768688175870E-11 
WEIGHT(3)=  1.758405638619854E-08 
WEIGHT(4)=  1.371239148976848E-06 
WEIGHT(5)=  4.435096659959217E-05 
WEIGHT(6)=  7 . 15550107S431907E-04 
WEIGHT(7)=  6.488956601264211E-03 
WEIGHT(8)=  3 . 644041585109798E-02 
WEIGHT(9)=  1.439979241604145E-01 
WEIGHT( 10)=  8 . 123114134235980E-01 
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SUM  =  0. 

ZETA  =  2 . * (Z*CDSQRT (Z) )/3 . 


100 


DO  100  1=1,  10 

SUM  =  SUM  +  (WEIGHT ( I) / ( 1 . 0+ (ZEROES (I) /ZETA) ) ) 
CONTINUE 

GQA  =  0 . 5*DSQRPI/ (CDSQRT (CDSQRT(Z) ) ) ♦SUM 

RETURN 

END 
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